Thermodynamics and the Global Optimization 
of Lennard-Jones clusters 
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Theoretical design of global optimization algorithms can profitably utilize recent statistical me- 
chanical treatments of potential energy surfaces (PES's). Here we analyze the basin-hopping al- 
gorithm to explain its success in locating the global minima of Lennard-Jones (LJ) clusters, even 
those such as LJ38 for which the PES has a multiple-funnel topography, where trapping in local 
minima with different morphologies is expected. We find that a key factor in overcoming trapping is 
the transformation applied to the PES which broadens the thermodynamic transitions. The global 
minimum then has a significant probability of occupation at temperatures where the free energy 
barriers between funnels are surmountable. 



I. INTRODUCTION 



Global optimization is a subject of intense current in- 
terest, both scientific and commercial. For example, im- 
proved solutions to optimization problems, such as the 
travelling-salesman problem and the routing of circuitry 
on a chip, could lead to cost reductions and improved 
performance. El In the chemistry and physics communi- 
ties motivation is often provided by the structural insight 
which can be derived from finding the lowest energy con- 
figuration of a (macro)molecular system. Therefore, the 
development of better global optimization algorithms is 
an important task. It is a challenge which needs to be 
guided by theoretical insight rather than proceeding on 
purely empirical or intuitive grounds. 

The main difficulty associated with global optimization 
is the exponential increase in the search space with sys- 
tem size. One expression of this problem is Levinthal's 
paradox which points out the apparent impossibility of 
finding the native state of a protein.0 The number of pos- 
sible conformations of a typical protein is so large that it 
would take longer than the age of the universe to find the 
native state if the conformations were sampled randomly, 
even at an unfeasibly rapid rate. This problem can be 
more rigorously defined using computational comnlexity 
theory£l Finding the global minimum of a proteinu or a 
clusterO is 'NP-hard', a class of problems for which there 
is no known algorithm that is guaranteed to find a solu- 
tion in polynomial time. 

However, we know that proteins are able to reach 
their native state relatively rapidly. The resolution of 
Levinthal's paradox lies in the fact that the search is not 
random — the paradox results from the implicit assump- 
tion that the potential energy surface (PES) is flat. In re- 
ality, the topography of the PES is crucial in deternwiing 
the ability of a system to reach the global minimum fl For 



example, the system is thermodynamically more likely 
(than random) to be in those regions of the PES that 
are low in energy, and downhill pathways can guide the 
system towards certain configurations. A funnel, a set of 
downhill pathways that converge on a single low energy 
minimum, combines these two effects. It has been sug- 
gested that the PES's of proteins are characterized by a 
single deep funnel and that this feature undedies the abil- 
ity of proteins to fold to their native state.Q-a Indeed it 
is easy to design model single- funnel PES's that result in 
efficient relaxation to the glahaL minimum, despite very 
large configurational spaces.tlEilj 

Therefore, global optimization should be relatively 
easy for those systems with a single-funnel PES topog- 
raphy. However, there are many systems for which the 
topography of the PES does not permit escape from the 
huge number of disordered configurations. These systems 
are more likely to form glassy or amorphous structures 
on cooling. In such cases it is very difficult for global 
optimization to succeed. This is especially true if the 
optimization method, e.g. simulated annealing, tries to 
mimic some physical process, because the natural time 
scales for the formation of order are much longer than 
those that can be probed by computer simulation. 

It has been suggested that the PES's of such 'glassy' 
materials are rough, with many funnels leading to dif- 
ferent low enei:gjy|-ajnorphous configurations, rather than 
to the crystal.liinlS However, one does not need to have 
much complexity to make global optimization difficult; 
two funnels are enough if the system is more likely to en- 
ter the secondary funnel on descending the PES.EI In this 
paper we examine one such example from the realm of 
atomic clusters and compare it with another much eas- 
ier problem. Our aims are to pinpoint why many global 
optimization algorithms arc likely to fail for this cluster, 
to understand the reasons for the success of the 'basin- 
hopping' method and to deduce lessons for the design of 
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global optimization algorithms. In particular, we high- 
light the influence of the thermodynamics of the clusters 
on these questions. A brief [description of some of the 
results has already appeared.Ej 



II. LENNARD-JONES CLUSTERS 

The Lennard- Jones (LJ) potential,Ei which provides a 
reasonable description of the interactions between rare 
gas atoms, is given by 
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(1) 



where e is the pair well depth and is the equilib- 

rium pair separation. LJ clusters have been much used 
as a test system for global optimization methods that 
are designed for configurational problems, and it is likely 
that nearly all the, global minima up to 150 atoms have 
now been found.EjEj Most of the global minima hajte 
structures that are based upon the Mackay icosahedra,EZI 
but there are a number of exceptions. The LJag global 
minimum is a face-centred cubic (fee) truncated octa- 
hedron (Fig Ha), and for 75-77 and 102-104 atopis the 
global minima are based on Marks decahedra.Ea These 
clusters are interesting because the global minima are 
much more difficult to find by an unbiased global op- 
timization method. They were all initially discovered 
by construction using physical insight Since then 
the truncated ■octahedron . has been found by a num- 
ber of methods J1j'E3"oEjcJ The Marks decahedral global 
minima have |Ouly been found by the 'basin-hopping' 



approach, lSEJ the method we analyse in this paper, and 
a modified genetic algorithm which searches the same 
transformed surface. l3 




For LJ38 the second lowest energy minimum is an in- 
complete Mackay icosahedron with point group sym- 
metry (Figure |l|b), which lies only 0.676e higher in en- 
ergy than the fee global minimum. As might be expected 
from their structural dissimilarity, the two lowest energy 
minima are well separated on the PES. Using a method 
which walks over the^EES from minimum to minimum 
via transition states,L3~L3 we have been able to find paths 
connecting these two structures. The Lowest energy path 
that we found is depicted in Figure |2[Ej It clearly shows 
that there is a large activation energy associated with 
passage between the two lowest energy minima. The fee 
and icosahedral regions of configuration space represent 
two distinct funnels on the PES. The minima at the top 
of the barrier are from the region of configuration space 
associated with the liquid-like state (Figure ^). At tem- 
peratures below the melting point the Boltzmann weights 
of these states are small implying that there is -also a 
large free energy barrier between the two funnels .E3 The 
dynamics confirm this picture; below the melting point 
a simulation run started in one of the funnels always re- 
mains trapped in that funnel. 
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FIG. 2. Pathway between the fee global minimum of 
LJ38 and the lowest energy icosahedral minimum. The tran- 
sition states on the pathway are marked by diamonds and 
the minima by circles. 



FIG. 1. (a) The LJsg global minimum is an fee truncated 
octahedron {E = -173.928427e). (b) Second lowest energy 
minimum of LJ38 {Ec — — 173.252378e). The structure is an 
incomplete Mackay icosahedron. (c) The LJ55 global mini- 
mum is a Mackay icosahedron {Ec = — 279.248470e). 

For the above reasons we choose to examine in 
more detail the behaviour of LJ38, and as a contrast 
LJ55. The tJss global minimum is a complete Mackay 
icosahedroncZl (Figure |l|c), a structure that is easily found 
by any reasonable global optimization algorithm. 



The shapes of the fee and icosahedral funnels are also 
rather different. For the fee funnel there is a large en- 
ergy gap (2.G72e) between the truncated octahedron and 
the next lowest energy minimum. In contrast there are 
many low energy icosahedral structures; indeed there are 
at least 24 icosahedral minima lower in energy than the 
second lowest energy fee minimum (Figure 0). The fee 
funnel is narrow, whereas the icosahedral funnel is wider 
and has a broad, fairly flat bottom. 



2 



0.9 




Energy / £ 

FIG. 3. Probability distribution for the potential energy 
of LJ38 minima obtained by systematic quenching from a mi- 
crocanonical molecular dynamics (MD) run at E = — 148.2e. 
866 distinct minima were obtained from the 2500 quenches. 
We have divided the minima up into three sets by their en- 
ergy, i.e. the Oh global minimum, icosahedral minima and 
minima that are associated with the liquid-like state. 

The LJ55 PES is very different. All the low energy 
minima for LJ55 are based on the Mackay icosahedron. 
The first three peaks above the global minimum in the 
probability distribution of minima in Figure ^ represent 
Mackay icosahedra_with one, two and three surface de- 
fects, respectivelyra Furthermore, the Mackay icosahe- 
dron is 2.644e lower than any other minimum. The PES 
of LJ55 has a single deep funnel. 
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FIG. 4. Probability distribution for the potential energy 
of LJ55 minima obtained by systematic quenching from a mi- 
crocanonical MD run at E = — 209e. 1153 distinct minima 
were obtained from the 1191 quenches. We have divided 
the minima up into five sets by their energy: the Mackay 
icosahedral global minimum, Mackay icosahedra with differ- 
ing number of defects (M-|-nd) and those that are associated 
with the liquid-like state. 



III. THERMODYNAMICS OF LJgs AND LJ55 

To calculate the thermodynamic propedies of our clus- 
ters we use the superposition methodpjO in this ap- 
proach the density of states, il{E), is constructed using 
information about a set of minima on the PES (such as 
those represented in Figures || and ||). This choice is 
motivated by two considerations.|-J&stly, methods such 
as multi-histogram Monte CarloJi^'eS often the method 
of choice for clusters, are not able to describe the low 
temperature thermodynamics for LJ38 correctly because 
in simulations the cluster is not able to cross the large 
free energy barrier between the fee and icosahedral fun- 
nels. It might be possible to overcome this difficulty bK 
using more advanced techniques such as jump- walkings 
or umbrella sampling,Cj however these methods can be 
computationally demanding. The advantage of the su- 
perposition method is that it can be used to calculate 
absolute densities of states for different regions of phase 
space, and therefore there is no need to cross any free 
energy barrier. Secondly, the superposition method can 
give greater physical insight into the relationship between 
the topography of the PES and the thermodynamics. 

We review the superposition method here because it 
will be important later when we formulate the thermo- 
dynamic properties of a transformed energy surface, and 
because some additional developments need to be made 
before it can be applied to LJag. The fundamental rela- 
tion is 

n{E) = J2 ^s^siE), (2) 

where the sum is over all the geometrically distinct min- 
ima on the PES, ^IsiE) is the density of states of a single 
minimum s, Eg is the potential energy of minimum s, 
and Us, the number of permutational isomers of mini- 
mum s, is given by Us = 2N\/hs where hs is the order 
of the point group of s. This equation is exact and just 
expresses the division of configuration space into basipa 
of attraction that surround each minimum on the PES.lj 
One difficulty with the above equation is that ^s{E) is 
not known a priori. However, if the basin is assumed 
to be harmonic then the form of Vls{E) is simple, and 
using this expression a qualitative picture of the thermo- 
dynamics can be obtainedEj Furthermore, anharmonic 
expressions for ^s{E) have also been derived which are 
able to describe the thepjnodynamics of Lennard- Jones 
clusters more accurately.c3 

The second difficulty with equation (||) is that for all 
but the very smallest clusters the sum involves an impra|C=. 
tically large number of minima. Hoare and Mclnnes,c£l 
and more recently Tsai and Jordan,EZI have enumerated 
lower bounds to the number of geometric isomers for LJ 
clusters from 6 to 13 atoms. This number rises exponen- 
tially with N . Extrapolating the trend gives for LJ55 an 
estimate of 10^^ geometric isomers. In such a case, as it 
is not possible to obtain a complete set of minima, one 
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instead has to use a representative sample, which, for ex- 
ample, could be obtained by performing minimizations 
from a large set of configurations generated by a molec- 
ular dynamics trajectory. However, one needs a way to 
correct for the incomplete nature of the sample. This cor- 
rection can be achieved by weighting the density of states 
for each known minimum by gs , the number of minima of 
energy ~ Eg for which the minimum s is representative. 
Hence, 



(3) 



E^<E 



where the prime indicates that the sum is now over an in- 
complete but representative sample of minima. The effect 
of gs can be incorporated by using ps{E'), the probabil- 
ity of obtaining s in a minimization from a configuration 
generated by a microcanonical simulation at energy E' . 
If the system is ergodic on the time scale of the simula- 
tion, then ps{E') oc gsnsft{E')s- Hence, 



(4) 



Since we know the low energy limiting form of r2(i?), 
where the contribution of the lowest energy minimum 
dominates, the proportionality constant in (4) can also 
be obtained. 

This rewei ghtigg .technique is analogous to histogram 
Monte CarlOjCilBcj but instead of determining the con- 
figurational density of states from the canonical potential 
energy distribution, g, effectively a density of minima, is 
found from the microcanonical probability distribution of 
being in the different basins of attraction. Hence there 
are some similarities with the microcanonical multihis- 
togram approach of Calvo and Labastie.L^I 

The accuracy of the method depends on the statistical 
accuracy oips{E')\ this probability distribution needs to 
be obtained at an energy where all relevant minima are 
significantly sampled. This is possible for LJ55;e2I how- 
ever, for LJ38 there is no energy at which both the fee 
and icosahedral minima can both be sampled because of 
the large free energy barrier. Therefore, in the latter case 
the density of states is obtained by adding two terms: 



n{E) = ni,,{E) + n,,,t{E). 



(5) 



The density of states of the fee region of phase space, 
^^fcc(^'), can be obtained by summing the density of 
states for all the low energy fee minima using equa- 
tion (^. (In fact the summation can be dispensed with 
since the term from the global minimum is the only term 
in Q.icc{E) that ever contributes significantly to Q.{E).) 
At the melting point an LJ38 cluster samples both the 
icosahedral and the liquid-like regions of configuration 
space and so equation (||) can be used to find forest 
Once ^1{E) has been determined, various thermodynamic 



quantities can be calculated by 
dard thermodynamical formulae.L 
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FIG. 5. Equilibrium thermodynamic properties of the 
untransformed LJ38 PES in the canonical ensemble, (a) 
The probability of the cluster being in the fee, icosahedral 
and 'liquid-like' regions of bound configuration space, (b) 
The heat capacity, Cv The sample of 661 minima was ob- 
tained from 2500 quenches along a microcanonical MD run 
a.t E = -150.1e. 



The resulting equilibrium thermodynamic properties 
of LJ38 are shown in Figure ^. The anharmonic 
expression |Jpr Sis was used with the anharmonicity 
parametersEj determined by the technique give in ref. 
As there are many low energy icosahedral minima and as 
these minima have a lower mean vibrational frequency 
than the truncated octahedron, the configurational en- 
tropy associated with the icosahedral funnel is larger than 
for the fee funnel. As a result of this entropy difference 
there is a transition from the fee to the icosahedral re- 
gions of configurational space at low temperature. The 
transition is centred at a temperature of ~ 0.11 efc~^ and 
gives rise to the small peak in the heat capacity (Figure 
^d) . The finite-size analogue of the bulk melting transi- 
tion occurs at ^ Q_17efc~^, producing the main peak in 
the heat capacity.o These transitions are also reflected 
in the occupation probabilities for the different 'phases' 
(Figure ||a) . 
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We can understand the effect these thermodynamic 
transitions have on global optimization by considering 
the result of cooling a cluster from the liquid-like state. 
The cluster is thermodynamically much more likely to 
enter the icosahedral funnel on quenching because the 
free energy associated with the icosahedral structures at 
the melting point is lower than that associated with the 
fee structures. Furthermore, this effect is reinforced by 
the topology of the PES. The structure of simple pOlxmiic 
liquids has significant polytetrahedral character;LJO in 
contrast fee structures have no polytetrahedral charac- 
ter, and icosahedra are somewhere in between. As a 
result of this greater structural similarity, the icosahe- 
dral funnel is more accessible from the liquid-like state. 
(A similar effect has been noted by Straley who found 
that crystallization from a simple liquid is much easier 
in a curved space, where a polytetrahedral packing is the 
global minimum, than in norpiaLflat space where a close- 
packed crystal is most stable.E3C3) The combined effect of 
these two features is to make it extremely probable that 
on cooling the cluster will enter the icosahedral funnel. 
On further cooling the cluster will remain trapped in this 
funnel even when the fee structures become lower in free 
energy because of the large free energy barrier between 
the two funnels. Clearly this behavior will present diffi- 
culties for annealing methods, which simply try to follow 
the free energy global minimum down to zero tempera- 
ture by gradual cooling. 

The only opportunity to find the truncated octahedron 
by conventional simulation is at high temperature in the 
narrow window where the fee structures have a small, 
but not yet insignificant, probability of being occupied 
and where the Boltzmann weights for the intermediate 
states, which mediate transitions between the fee and 
icosahedral funnels, are large enough to make the free 
energy barrier surmountable. Indeed, we did observe the 
truncated octahedron in the microcanonical simulation 
(T = 0.185eA:^^) from which we obtained the distribution 
of minima depicted in Figure ^. However, minimization, 
only led to this structure once in the entire 0.25 fis run.Lj 

Some thermodynamic properties of LJ55 are shown in 
Figure |^. The partition function for this system was ob- 
tained previcmsly in the development of the superposi- 
tion methodjfj and it gives very good agreement with 
thermo dyna mic properties derived by more conventional 
means.E^S There is a single peak in the heat capacity 
curve corresponding to the melting transition at which 
the cluster passes from the Mackay icosahedron (per- 
haps with one or two defects) into the liquid-like state. 
For simulations in the middle of this melting region 
the cluster passes back and forth between the different 
states.E3c3 Similarly, on cooling from the liquid the clus- 
ter enters the icosahedral funnel relatively easily. LJss's 
single funnel PES presents no thermodynamic obstacles 
to global optimization. 
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FIG. 6. Equilibrium thermodynamic properties of the 
untransformed LJ55 PES in the canonical ensemble, (a) The 
probability of the cluster being in the Mackay icosahedron, 
with one or two defects, and 'liquid-like' regions of bound 
configuration space, (b) The heat capacity, C„. 



IV. GLOBAL OPTIMIZATION BY 
BASIN-HOPPING 

The global optimization method that we analyze here 
belongs jta the 'hypersurface deformation' family of 
methods.Eil In this approach the potential energy func- 
tion is transformed in a way that is hoped to make 
global optimization easier. The transformations usually 
try to lower the number of local minima — thus reduc- 
ing the search space — and/or the barrier heights between 
minima — thus increasing the rates of transitions between 
minima. However, little attention is usually paid to the 
effects of the transformation on the thermodynamics of 
the system. 

Once the global minimum of the transformed PES is 
found it is mapped back to the original surface in the 
hope that this will lead to the global minimum on the 
original PES. However, there is no guarantee that the 
global minima on the two surfaces are related and of- 
ten there are good reasons to think that they are not. 
For example, one method suggested-for clusters is to in- 
crease the range of the potentialHSS Such a transforma- 
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tion can dramatically reduce the rLtujnbfir of minimacJ 
and also lower the barrier heights Eji£3 However, it has 
been shown that changing the range of the potential for 
clusters can leacL-tt] Bjia| ny changes in the identity of the 
global ininimum.E3'E!}E3 

In contrast, the transformation that we apply to the 
PES is guaranteed to preserve the identity of the global 
minimum. The transformed potential energy Ec is de- 
fined by: 



(6) 



where X represents the vector of nuclear coordinates and 
min signifies that an energy minimization is performed 
starting from X. Hence the energy at any point in con- 
figuration space is assigned to that of the local minimum 
obtained by the minimization, and the transformed PES 
consists of a set of plateaus or steps each corresponding 
to the basin of attraction surrounding a minimum on the 
original PES. 




FIG. 7. A schematic diagram illustrating the effects of 
our potential energy transformation for a one-dimensional 
example. The solid line is the potential energy of the orig- 
inal surface and the dashed line is the transformed energy. 

One can easily see that the transformation will have a 
significant effect on the dynamics. Transitions to a lower 
energy step are barrierless, and so relaxation down a fun- 
nel is like descending a multi-dimensional staircase (Fig- 
ure Furthermore transitions can occur at any point 
along the boundary between basins of attraction, whereas 
transitions on the untransformed surface can occur only 
when the system passes along the transition state valley. 
As a result intrawell vibrational motion is removed and 
the system can hop directly between basins of attraction 
at each simulation step, hence our name for the algo- 
rithm, 'basin-hopping'. In structural terms the system 
can now pass through large repulsive barriers on the un- 
transformed PES. Similarly, if the method were applied 
to a polymer system it should remove entanglement ef- 
fects because chains can pass through each other. 

However, the transformation does not entirely remove 
the energy barriers between two funnels but only the 
component due to transition states (Figure 0). There- 
fore, it is not self-evident that it will aid global opti- 
mization on multiple-funnel PES's; this effect depends 



on how the transformation affects the thermodynamics 
of the system. 

To explore the transformed PES we simply use a 
canonical Monte Carlo simulation. For clusters we also 
have to consider how to restrict our search of configu- 
ration space to that for bound clusters. This problem 
is more pressing than for the original PES, since the 
transformation has removed many of the energy barri- 
ers to the dissociation of the cluster. We have considered 
two approaches to achieve this. The first is to place the 
cluster in a tight-fitting box; the second is to reset the 
coordinates at each successful step to the relevant mini- 
mum. In the latter case the method is equivalent to Li 
and Scheraga's Monte Carlo Minimization approach,L3 
although they did not conceive of the method in terms 
of a transformation of the PES. White and Mayne have 
now shown that resetting the coordinates is probably the 
best strategyJI3 and it is the one we concentrate on here. 

We have found that the 'basin-hopping' algorithm per- 
forms well. In our application to Lennard-Jones clus- 
ters it found all the known lowest energy configurations 
up to 110 atoms, including those with non-icosahedral 
global minima.E^ Furthermore, it has performed impres- 
sively for a wide range of cluster systems some of which 
have a considp¥|arh]AL| more rugged PES than Lennard- 
Jones clusters. E30'E3 

In Figure |^ we show some basin-hopping trajectories 
for LJ38 and LJ55. It can be seen that the cluster is able 
to pass between the fee and icosahedral funnels of LJ38 
over a range of temperature. Therefore, the free energy 
barriers must be lowered by the transformation. In the 
next section we will examine the thermodynamics of the 
transformed PES to discover the basis for this change. 
The predicted acceleration of the dynamics is also evi- 
dent from the frequency of transitions between the differ- 
ent states of LJ38 and LJ55 (Figure^). The timescale (in 
numbers of steps) for these processes is much more rapid 
than for similar processes on the original PES, where 
even diffusion across a flat free energy landscape can be 
slow. However, it should be remembered that a single 
MC step on Ec is much more computationally expensive 
than an MD step or an MC cycle on the untransformed 
PES, because it involves a minimization. 



A. Thermodynamics for Ec 



To calculate the thermodynamics for the transformed 
energy surface, E^ the partition of configuration space 
that we used in section HI (equation (|^)) is appropriate 
since the plateaus of Ec correspond to the basins of at- 
traction on the original PES. However, we need to derive 
the form of flgiE). The configurational density of states 
of a plateau s, flcsiEc) is simply 



nc,siEc) = 6iEc-E,)As 



(7) 



6 



where Ag is the hyperarea of the basin of attraction of 
minimum s. This expression is then convoluted with the 
kinetic density of states to give 

where m is the mass of an atom, h is Planck's constant, 
K, the number of vibrational degrees of freedom, is iN — Q 
and 9 is the Heaviside step function. As As cannot be 
readily calculated, we cannot use equation (H) to obtain 
the total density of states but instead employ equation 
(^. For Ec the latter approach can even be applied to 
LJ38 since, as Figure ||b shows, conditions can be found 
where the fee, icosahedral and liquid-like regions of con- 
figurational space are all significantly sampled. How- 
ever, since we do our simulations on Ec using standard 
Metropolis Monte Carlo sampling, we actually use the 
canonical analogue of equation (^) : 

Z(/3)(x^p,(/3')J^, (9) 

where Z is the partition function, the inverse tempera- 
ture (3 ~ 1/kT and the probability distribution was ob- 
tained from a basin- hopping trajectory at /?'. To use 
this equation we must first Laplace transform the ex- 
pression for ils in (8) to give Zg, the partition function 
of a plateau: 

Equation (||) then reduces to 

Z(/3)(x^p,(/3')e-^=(''-^'). (11) 

s 

It is interesting to note that t hi s t^ qu ation is identical 
to that used in histogram MC.cilHu^ This equivalence 
stems from the fact that the potential energy is constant 
on each plateau s. Hence the multi-histogram techniques, 
in which information from a number of runs at different 
temperatures is used to construct the partition function, 
could also be applied. However, such an approach was 
not necessary in this study. 

The As values depend upon how configuration space 
is restricted to bound clusters. When resetting the co- 
ordinates to those of the minimum at each successful 
step, configuration space is restricted to a hypersphere 
(with radius corresponding to the maximum step size) 
around each local minimum, and so the thermodynanues 
(and the performance of the basin-hopping algorithmO) 
is somewhat dependent on the maximum step size used. 
When using a container the total accessible configuration 
space, '^HsAs, is the volume of a K-dimensional hyper- 
sphere whose radius corresponds to that of the container. 
Therefore, the thermodynamic contribution of irrelevant 



regions of configuration space corresponding to clusters 
of low density, or with atoms evaporated, increases with 
container size. Hence, it is advantageous to use a tight- 
fitting container, and for this situation the thermody- 
namics are similar to when the coordinates are reset. 
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FIG. 8. Ec as a function of the number of steps in 
basin-hopping runs for LJ38 at (a) T = 0.4efc"^ and (b) 
T = 0.9ek~^ In (c) we show a close-up of one section of (b) 
in which the system passes from the icosahedral region of 
configuration space to the fee region and then back again, 
(d) A basin- hopping run at T = lO.Oefe"^ for LJ55. 

One further consideration is that one can imagine sit- 
uations where resetting the coordinates breaks detailed 
balance. For two minima A and B which are adjacent on 
the PES, it may be that no part of the basin of attraction 
of minimum B is within the maximum step size of mini- 
mum A, whereas the basin of attraction of A lies within 
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the maximum step size of minimum B. If so the system 
can pass from B to A but never from A to B, thus break- 
ing the detailed balance condition. Although this means 
that the 'resetting' implementation of the basin-hopping 
algorithm might not formally produce a canonical en- 
semble, in practice it does to a good approximation — 
the thermodynamic properties calculated using equation 
(|ll|) for samples generated at different temperatures are 
in reasonable agreement. 




temperature / e k"' 

FIG. 9. Equilibrium thermodynamic properties of the 
transformed LJss PES when the maximum step size is fixed 
at 0.35a. (a) The probability of the cluster being in the fee, 
icosahedral and 'liquid-like' regions of bound configuration 
space, (b) The configurational component of the heat ca- 
pacity. These properties were Calculated from a probability 
distribution produced by a 500000 step run at T = 0.9efc~^. 

For the transformation of the LJ38 PES to reduce the 
free energy barriers between the funnels, the probability 
of the system occupying the intermediate states between 
them must be non-negligible under conditions where the 
icosahedral and fee structures also have significant occu- 
pation probabilities. The thermodynamics of theplxans- 
formed PES have just these properties (Figure The 
transitions have been smeared out and there is now only 
one broad peak in the heat capacity. Most significantly, 
the probability that the system is in the basin of attrac- 
tion of the global minimum decays much more slowly, 



and the temperature range for which both the liquid-like 
minima and the global minimum have significant prob- 
abilities is large. Basin-hopping runs anywhere in this 
temperature range are able to locate the global mini- 
mum rapidly from a random starting point. At the lower 
temperatures in this range, e.g. T — 0.4efc^^, there is 
still a small free energy barrier and so the cluster can be 
trapped in one of the funnels for many steps (Figure ^), 
however at higher temperatures the system passes more 
rapidly between the states (Figure ||b and c). 



g_ 8000- 
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0.4 0.6 0.8 1 1.2 L4 1.6 1.8 2 

temperature /£k"' 

FIG. 10. Mean first-passage time for the basin-hopping 
algorithm to reach the LJ38 global minimum as a function of 
temperature when the maximum step size is fixed at 0.35(j. 
Each point represents an average over 100 runs starting from 
a random configuration. 

In this way we can correlate the performance of the 
basin-hopping algorithm with the thermodynamics. The 
temperature dependence of the first-passage time for 
reaching the global minimum from a random configu- 
ration is shown in Figure As the temperature de- 
creases the probability of the system being in the high 
energy intermediate states between the two systems de- 
creases. The resulting increasing free energy barrier 
causes the first-passage time to rise at low tempera- 
tures. Interestingly, and perhaps unexpectedly, there is 
little rise-at-Jriigher temperature. For many systems, e.g. 
proteins the first-passage time has a minimum as a 
function of temperature; the rise at high temperature 
is because the equilibrium probability of being in the 
global minimum tends to zero. However, for the basin- 
hopping algorithm the first-passage time at high temper- 
ature is approximately constant E3 In fact the probability 
of being in the global minimum never goes to zero, e.g. 
po^{T = 00) = 0.0024. 

The thermodynamics on the transformed surface for 
LJ55 are also considerably broadened and the transitions 
have been smeared out over a large temperature range 
(Figure [Tl| ). Consistent with this, there is a greater 
probability that one of the states intermediate between 
the global minimum and the liquid is occupied. As 
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expected for the single funnel topography of LJ55, the 
basin-hopping algorithm needs far fewer steps to find the 
global minimum than for LJsg. It is remarkable that for 
a wide range of temperature (T > 1.5efc~^) the method 
on average requires fewer than 200 minimizations to find 
the global minimum from a random starting point (Fig- 
ure |2|) despite the estimated 10^^ minima on the PES. 
This is a testament not just to the efficiency of the basin- 
hopping algorithm but also to the dramatic effect that 
a funnel has in guiding the system towards the global 
minimum. As with LJ38 the first-passage time is ap- 
proximately constant at high temperature (Figure 12); 

P7,(r = cx)) = 0.11. 



Us exp{-PEs 



(12) 




M+3d 




temperature /£k 

FIG. 11. Equilibrium thermodynamic properties of the 
transformed LJ55 PES when the maximum step size is fixed 
at 0.35a. (a) The probability of the cluster being in the 
Mackay icosahedron, with one to three defects, and in 'liq- 
uid-like' regions of bound configuration space, (b) The con- 
figurational component of the heat capacity. The properties 
were calculated from a probability distribution produced by 
a 100000 step run at T = lO.Oefc"^ 



The superposition method allows us to connect the 
thermodynamics to the properties of the PES. Hence we 
can understand the thermodynamics of the transformed 
surface by examining expressions for the canonical prob- 
ability that the system is in a minimum or on a plateau 
s. On the untransformed PES 



where Vs is the geometric mean vibrational frequency 
of minimum s, and we have used the harnmnic 
approximation — the more accurate anharmonic formCj is 
complicated and does not provide any additional physical 
insight. For the transformed PES 



UsAs exp(-/3_Eg) 



(13) 



In equations ( p^ ) and (^3|) the sums are over all the min- 
ima on the potential energy surface. 



5 10 15 20 

temperature /£k"' 

FIG. 12. Temperature dependence of the mean 

first-passage time for the basin-hopping algorithm to reach 
the LJ55 global minimum when the maximum step size is 
fixed at 0.35cr. Each point represents an average over 100 
runs starting from a random configuration. 



The two expressions (12) and ( p^ ) differ only in the vi- 
brational frequency and As terms. The fee to icosahedral 
and the icosahedral to liquid transitions are caused by 
the greater number of minima (and therefore the larger 
entropy) associated with the higher temperature states. 
On the untransformed surface this effect is reinforced by 
the decrease in the mean vibrational frequencies with in- 
creasing potential energy (Figure |l3|a) . Although the de- 
pendence of Vs on Es may seem relatively small, because 
Vs is raised to the power 3A^ — 6 it leads to a signifi- 
cant increase in the entropy of the higher energy states, 
sharpening the transitions and lowering the temperature 
at which they occur. In contrast, Ag decreases rapidly 
with increasing potential energy, (Figure p^). The de- 
crease in As reduces the entropy of the higher energy 
states, causing the transitions to be broadened and the 
temperature at which they occur to increase. 
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FIG. 13. (a) Geometric mean vibrational frequencies, Vs, 
and (b) basin of attraction areas, As, for a sample of LJss 
minima. 



To obtain the Ag values shown in Figure [T^ b we 
used Ps{l3'), the probabihty with which a plateau s 
was visited in a basin-hopping trajectory. As PsiP) cx: 
gsUgAs exp{—(3Es), relative values of gsAg can be found. 
A similar inversion of ps{E') for the |Untransfornied sur- 
face allows values of gs to be found. c2l However, the gs 
and gsAs values are for different sets of minima. The two 
samples of minima overlap completely for the low energy 
minima, and so A^ values for these minima can be cal- 
culated directly. For higher potential energies, where the 
spectrum of minima is quasi-continuous, we find the av- 
erage value of As in an interval of the potential energy, 
AEc, from gsAg/ 9s where the sums are only over 
the minima in the two samples which occur in that energy 
interval. 

The superposition approach also allows us to com- 
ment on the thermodynamics of a related optimization 
algorithm which steps directly between connected min- 
ima on the PES (the steps arepachieved using a transi- 
tion state searching algorithm) E3~E3 As with the basin- 
hopping algorithm the detrimental effects of vibrational 
motion are thus removed. This method has been suc- 



cessfully used to fincl low, energy structures of amor- 
phous semiconductorsOL^iIa and its performance fps a 
few Lennard- Jones clusters has also been illustrated. 13 If 
each step is accepted with a Metropolis criterion then it 
is easy to formulate the thermodynamics for this discrete 
space of minimajlj 



Us exp(-/3gs) 



(14) 



This expression is intermediate between equations ( p^ ) 
and (|l3). The thermodynamics have been broadened 
with respect to the original PES, but not as much as 
for the 'basin- hopping' transformation of the PES. The 
above approach is probably less efficient than 'basin- 
hopping' because of the expense of searching for a tran- 
sition state at each step. However, it has the advantage 
that interesting information about the. PES, .yuph as re- 
action pathways for complex processesHEIEj'Ej and the 
connectivity patterns between the minima,L3 can also be 
obtained. 



V. CONCLUSIONS 

We can now explain in detail why the basin-hopping or 
Monte Carlo MinimizationEa method is successful. First, 
the staircase transformation removes the transition state 
barriers between minima on the PES without changing 
the identity of the global minimum, accelerating the dy- 
namics. Second, it changes the thermodynamics of the 
surface, broadening the transitions so that even for a 
multiple- funnel surface such as that of LJ38, the global 
minimum has a significant probability of occupation at 
temperatures where the free energy barrier for passage 
between the funnels is surmountable. 

It is this latter feature which is especially important in 
overcoming the difficulties associated with multiple fun- 
nels and represents a new criterion for designing success- 
ful global optimization methods. Most previous hyper- 
surface deformation schemes have been developed with- 
out regard to the thermodynamic effects of the transfor- 
mation and so, in some cases, they may make optimiza- 
tion more difficult. The use of Tsallis statisticsEjEZl to 
improve annealing algorithms is another example of how 
thermodynamic insight can be used. 

The broadened transition also means that the global 
minimum can be found by simulations over a relatively 
broad range of temperature. Therefore, simple canoni- 
cal Monte Carlo is sufficient to search the transformed 
surface and the performance is relatively insensitive to 
the choice of temperature. However, any of the more 
advanced techniques for searching PES's such as simjtu 
lated annealing,EJ_ffenetic algorithms,La jumn_*ralking,E2l 
Tsallis statistics,L3 and simulated temperingj£3 could be 
used to search the transformed PES. As many of these 
methods have been designed to work on untransformed 
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PES's it is not obvious that they will improve perfor- 
mance, and preliminary investigations have not yet in- 
dicated that any of them provides a significant advan- 
tage. It is also interesting to note that the most success- 
ful applications of genetic algorithms to clusters refine 
each new configurat inn generated by mutation or mat- 
ing by minimization£a'E3 this minimization step has been 
shown to be crucialo These genetic algorithms are there- 
fore searching the transformed PES, Ec, and it is prob- 
ably this feature which is responsible for their success. 
In fact, we have succeeded in obtaining comparable re- 
sults for Lpanard- Jones clusters with a modified genetic 
algorithm.L3 The most efficient method is likely to be 
system-dependent and draw upon features from a num- 
ber of different approaches. 

LJ75 could prove to be an interesting example to ex- 
plore the use of other techniques to search the trans- 
formed surface because its global minimum, even with 
the basin-hopping algorithm, is difficult to find. The free 
energy barriers between the decahedral and fee funnels, 
although reduced by the transformation of the PES, are 
still large enough to make global optimization difficult. 
Therefore, LJ75 may provide a suitable system in which 
to investigate algorithms that enhance barrier crossing 
or simulate an ensemble which has broader probability 
distributions than the canonical ensemble. 

The broadened transition that results from the stair- 
case transformation also differs markedly from the op- 
timum conditions for protein folding, if wc assume that 
proteins have evolved single-funnel surfaces in order to 
fold efficiently. A steep funnel provides a large ther- 
modynanHft, driving force for relaxation to the global 
minimuniQ'Q, and also leads to a sharp thermodynamic 
transition. However, in global optimizations one can- 
not make assumptions about the topography of the PES, 
and on a multiple-funnel PES features such as steepness 
can exacerbate the djfRculties associated with trapping 
in secondary funnelstl. 

Moreover, in protein folding there is the extra require- 
ment that the folded protein must remain localized in 
the native state, and this necessitates a sharp transition. 
There is no need for a global optimization method to 
mimic this property. Indeed, when applied to a PES 
with multiple funnels the optimum temperature for the 
basin-hopping approach is above the transition, where 
the system only spends a minority of its time associated 
with the global potential energy minimum. 
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